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Diffusion in coulomb crystals can be important for the structure of neutron star crusts. We 
determine diffusion constants D from molecular dynamics simulations. We find that D for coulomb 
crystals with relatively soft-core 1/r interactions may be larger than D for Lennard- Jones or other 
solids with harder-core interactions. Diffusion, for simulations of nearly perfect body-centered-cubic 
lattices, involves the exchange of ions in ring-like configurations. Here ions "hop" in unison without 
the formation of long lived vacancies. Diffusion, for imperfect crystals, involves the motion of defects. 
Finally, we find that diffusion, for an amorphous system rapidly quenched from coulomb parameter 
r = 175 to coulomb parameters up to T — 1750, is fast enough so that the system starts to crystalize 
during long simulation runs. These results strongly suggest that coulomb solids in cold white dwarf 
stars, and the crust of neutron stars, will be crystalline and not amorphous. 

PACS numbers: 66.30.-h 97.60.Jd 52.27.Lw 



I. INTRODUCTION 

Diffusion in coulomb plasma liquids has been well stud- 
ied pQ and is important for sedimentation of impurities 
in white dwarf (WD) gH3] and neutron stars (NS) [HI E] • 
Here ions, with a larger than average mass to charge 
ratio, sink in a strong gravitational field. This releases 
gravitational energy that can delay the cooling of metal 
rich WD [7J. However, we are not aware of numerical 
results for diffusion constants of coulomb crystals under 
Astrophysical conditions. Often the diffusion constant 
is simply assumed to be zero. This diffusion could be 
important for sedimentation in solid WD interiors, over 
long time scales, and for the structure of NS crusts. 

Solid diffusion can depend dramatically on the form of 
the interaction between particles and may be very slow 
for hard-core systems. For example, the binary Lennard 
Jones (LJ) system with a hard-core oc r~ 12 interaction 
forms a glass because of very slow diffusion [8]. In con- 
trast, the coulomb plasma with a soft 1/r core should 
have much faster diffusion. Therefore the Coulomb crys- 
tal may provide an important model system where diffu- 
sion is fast enough to be more easily studied by molecular 
dynamics (MD) simulations. 

In the laboratory, one can observe diffusion in complex 
(or dusty) plasma crystals. Complex plasmas (CP) are 
low temperature plasmas containing charged micropar- 
ticlcs, for a review see Fortov et al. [9j. Often the mi- 
croparticles are micron sized spheres that acquire large 
electric charges and the strong coulomb interactions be- 
tween microparticles can lead to crystallization. Indeed 
plasma crystals were first observed in the laboratory in 



1994 [TU]. Complex plasmas typically differ from White 
Dwarf interiors and Neutron Star crusts in a number of 
ways. First the microparticles feel additional fluctuating 
and friction forces because of interactions with the back- 
ground gas. Note that in stars, electron-ion interactions 
are small because of the large electron degeneracy. Sec- 
ond, the Dcbyc screening length A, sec Eq. [T] below, is 
often smaller in the CP than in a star (when measured 
in units of the lattice spacing). This changes the lat- 
tice type from body-centered-cubic (bcc) as expected in 
stars, to face-centered-cubic (fee) or other types in a CP. 
Finally in a CP there is an overall confining potential, 
and because of gravitational gradients it is often easier 
to study two-dimensional CP crystals. 

In two dimensions, one can have liquid, crystalline, 
and semi-crystalline states. Anomalous diffusion in 
semi- crystalline CP states has been observed at inter- 
mediate times [TT1 [T2"] . In anomalous diffusion the 
square of the displacement does not grow linearly with 
time. Langevin-dynamics simulations |13j find that 
microparticle-background gas interactions are important 
for this diffusion. 

The melting of colloidal crystal films has recently been 
studied [13] • Thick films (> 4 layers) were observed to 
melt at grain boundaries, while films 2 to 4 layers thick 
melted from both grain boundaries and from within crys- 
talline domains. We study diffusion at grain boundaries 
in Sec. iHlBl 

Three-dimensional CP crystals have been formed on- 
board the International Space Station under micrograv- 
ity conditions. Details of the experiment are presented 
in ref. |15j . The structural properties of the crystal were 
analyzed with bond angle metrics q^ and gg, see Section 
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IIIC Microparticles were found in regions with fee and 
hexagonal-close packing (hep) order (T6l [17], in agree- 
ment with MD simulations [15] . Khrapak et al. [18] 
studied freezing and melting of these CP crystals and 
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found diffusion to be relatively fast so that the system 
remained in equilibrium. Melting criteria for CP systems 
were presented by Klumov [TP] . 

We now focus on simple plasmas in three dimensions. 
The diffusion mechanism is interesting. Astrophysical 
systems are under great pressure that suppresses the for- 
mation of vacancies. Therefore diffusion, in a nearly per- 
fect crystal, should involve the exchange of neighboring 
ions. These exchanges, while common in some quantum 
systems, may be less common in classical systems. More 
complicated coulomb solids can involve a variety of dislo- 
cations, grain boundaries, and other imperfections. Dif- 
fusion in these systems probably involves motion of the 
imperfections, since this may be faster than particle ex- 
changes. Determining the diffusion constant for a system 
may help characterize the kinds and numbers of imper- 
fections. Note that the coulomb plasma has especially 
simple interactions. Therefore, it may be a very useful 
model system to study diffusion in the presence of com- 
plex imperfections. 

We emphasize that the coulomb plasma has no hard 
core interaction between ions, but only a relatively weak 
1/r repulsion. Therefore, it may be possible for ions to 
come relatively close to one another, if necessary for the 
motion of defects. This may be different from conven- 
tional condensed matter with hard cores. For example, 
MD simulations of defect motion in Magnesium focused 
on paths that involved only very small displacements of 
Magnesium atoms |20j . Imperfections may move much 
faster in a coulomb plasma. 

The motion of imperfections is important for equili- 
bration. For example, a coulomb liquid may freeze into 
an imperfect crystal state involving an excess of defects. 
There has been some work on nucleation in coulomb plas- 
mas, see for example [ST]. However present MD simula- 
tions of nucleation may have limitations from important 
finite size effects [22]. In this paper, we also study diffu- 
sion in amorphous systems to see if it is fast enough to 
allow crystallization. 

We focus on one component plasmas (OCP). We plan 
to study diffusion in multicomponent plasmas (MCP) in 
the future. As we discuss below, this may help address 
an important unsolved problem, the structure of MCP 
crystals. This is important for the thermal and electrical 
conductivity of NS crust [33] . Indeed X-ray observations 
of rapid NS crust cooling, after extended periods of accre- 
tion, strongly favor the formation of a crystalline rather 
than amorphous crust and may set limits on impurities 
[24H27j . In addition, pycnonuclear reactions, which are 
driven by quantum zero point motion at high densities, 
are exponentially sensitive to the structure of MCP crys- 
tals and the spatial locations of reactants [28]. These 
reactions may provide an important heat source in the 
crust of accreting NS [55]. Finally the distribution of 
dislocations, grain boundaries, impurities, and other im- 
perfections are important for mechanical properties of 
NS crust such as its breaking strain [30, 31 . The break- 
ing strain helps determine the maximum sized mountains 



that are possible on a NS, which are important for grav- 
itational wave radiation [30l [32]. The breaking strain 
also determines the maximum sized "star quake" that is 
possible. Sudden changes, or glitches, in the rotational 
period of pulsars [33j may involve crust breaking that 
could trigger the motion of superfluid vortices. In ad- 
dition Magnetar giant flares, extremely intense gamma 
ray flares from very strongly magnetized NS [34), may 
involve the catastrophic breaking of the crust because of 
very large magnetic stresses [35] . 

In previous work we determined liquid-solid phase 
equilibrium for a MCP system involving many ion species 
[55] , see also [37] ■ We performed a large scale MD sim- 
ulation where both liquid and solid phases were present. 
The solid phase in this simulation may have had a num- 
ber of imperfections. A knowledge of diffusion constants 
D may help determine the simulation time necessary for 
these imperfections to come into equilibrium. 

There have been previous calculations of D for coulomb 
liquids, starting with the MD simulations of Hansen et 
al. for the one component plasma (OCP) [l s . The one 
component plasma consists of ions, with pure coulomb 
interactions, and an inert neutralizing background charge 
density. Diffusion in the OCP in a strong magnetic field 
was considered by Bernu [38]. Hansen et al. have also 
calculated diffusion for binary mixtures [39 ■ 

Diffusion for a Yukawa fluid has been simulated by 
Robbins et al. [3D] and Ohta et al. [41] . In a Yukawa 
fluid ions interact via a screened coulomb potential My (r) , 

% (r) = ^e-^, (1) 

for two ions with charges Z\ and Zj, that are separated 
by a distance r. The OCP is equivalent to a Yukawa 
fluid, where all of the ions have the same charge and the 
screening length A is very large. 

The motion of ions in a WD or NS is largely classical 
because of their large mass. However at great densities, 
there could be quantum corrections that might increase 
D. These have been estimated for a liquid by Daligault 
and Murillo [32], and found to be very small. 

In this paper, we present classical MD simulations of 
one component crystals with Yukawa interactions in or- 
der to determine diffusion coefficients D. In Section ITT1 we 
describe our MD formalism and present results for diffu- 
sion coefficients in Section IIII1 We conclude in Section 
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II. FORMALISM 

We describe our MD simulation formalism. This is 
similar to what we used earlier to calculate D for liquid 
mixtures of carbon, oxygen, and neon [4] . We consider a 
one component system of oxygen ions where the ions are 
assumed to interact via screened Yukawa interactions, see 
Eq. [I] The Thomas Fermi screening length A, for cold 
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relativistic electrons, is 

A" 1 = 2a 1 ' 2 k F /K 1 ' 2 



(2) 



where the electron Fermi momentum kp is kp = 
(37r 2 n e ) 1 ^ 3 and a is the fine structure constant. The 
electron density n e is equal to the ion charge density, 
n e = Zn, where n is the ion density and Z is the 
ion charge. Our simulations are classical and we have 
neglected the electron mass (extreme relativistic limit). 
This is to be consistent with our previous work on neu- 
tron stars. However, the electron mass is important at 
lower densities in WD and this will decrease A. For rela- 
tivistic electrons, the ratio of A to the ion sphere radius 
a, 
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1/3 



(3) 



depends only on the charge Z and is independent of den- 
sity. For nonrelativistic electrons A/a can be somewhat 
smaller. In Section |III| we perform simulations for two 
values of A/a. 

The simulations can be characterized by a coulomb 
parameter T, 



2 D 2 



r = 



Z 2 e 



(4) 



Here T is the temperature. The system freezes near T = 
175 [43]. Note that this value of T may depend slightly 

on a mm]. 

Time can be measured in our simulations in units of 
one over the plasma frequency uj, p . Long wavelength fluc- 
tuations in the charge density can undergo oscillations at 
the plasma frequency. This depends on the ion charge Z 
and mass M, 



[4ne 2 Z 2 n 



M 



1/2 



(5) 



The diffusion constant D can be calculated from the 
velocity autocorrelation function Z(t), 



Z(t) 



i(t +t) -Vj-fto)) 
Vj(to) ■ Vj(i )) 



(6) 



where the average is over all ions j and over initial times 
<o- The velocity of the jth ion at time t is Vj(t). The 
diffusion constant is the time integral of Z(t), 



D = 



T 
M 



dtZ(t). 



(7) 



This Eq. works well to calculate D for liquids. However 
for crystals, D is smaller and the integral in Eq. [7] in- 
volves sensitive cancelations between regions where Z(t) 
is positive and negative. This makes Eq. [7] very difficult 
to use. 

Instead D can also be calculated from 



D(t) 



(|r 3 (t + < )~r,(to)| 2 ) 
6t 



(8) 



where the diffusion constant D is the large time limit of 
D(t), 



D = lim, 



(9) 



Here Tj (t) is the position of the jth ion at time t and the 
average in Eq. [8] is over all ions j and initial times £q- In 



principle, Eqs. 8[9 will have errors at large times t from 



the effects of periodic boundary conditions as \rj(t+to) — 
rj(t )\ becomes comparable to the size of the simulation 
volume. However diffusion is relatively slow so this is 
often not a problem until very large t. 

Note that D(t) can differ significantly from D for small 
t. For example, an ion undergoing thermal oscillations 
about an equilibrium lattice site will have Tj(t) — Tj(0) 
nonzero even though the ion remains near its original 
lattice site and there is no net contribution to diffusion. 
Therefore we define an alternative quantity D'{t) that 
has no contribution from ions that remain near their orig- 
inal lattice site, 



D\t) = 



6t 



(10) 

with t' = t + t . The cutoff radius R c is of order the 
lattice spacing, and will be discussed in Section |III| In 
the limit of very large times all ions move significantly so 
that D'{t) ->■ D(t) ast^oo. We observe that D'(t) is 
approximately independent of t, even for moderate t, so 
that 



D w D'(t) . 



(11) 



We use this equation, at finite t, to calculate D in Section 

IE 

The initial conditions are very important for determin- 
ing D because the system may contain different distribu- 
tions of defects and these distributions may take a very 
long time to equilibrate. We consider three classes of ini- 
tial conditions. The first class we call bcc and starts the 
ions with positions on a perfect body centered cubic (bcc) 
lattice and random thermal velocities. This may underes- 
timate the role of defects if there is not enough simulation 
time for thermal excitations to introduce an equilibrium 
distribution of defects. The second class of initial condi- 
tions we call imperfect crystal and starts the system from 
a liquid configuration that is cooled by rescaling the ve- 
locities until the system freezes. This may over estimate 
the role of defects if the system freezes into a very im- 
perfect state with more defects than would be present in 
thermal equilibrium. Note that imperfect crystal initial 
conditions may contain two or more micro-crystals with 
different orientations. Finally, we consider amorphous 
initial conditions where a liquid configuration is rapidly 
quenched to a much lower temperature. 

We evolve the system in time using the simple velocity 
Verlet algorithm [45]. We approximately maintain the 
system at constant temperature by simply rescaling the 
velocities every ten time steps. In Section [III| we present 
results for D. 
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III. RESULTS 



A. Body centered cubic lattice initial conditions 



We now present results for our MD simulations. We 
begin with a few results insensitive to initial conditions 
and then we discuss simulations with perfect lattice ini- 
tial conditions in Se ction |III A[ imperfect crystal initial 
conditions in Section IIIB[ and amorphous initial condi- 



tions in Section [ill C[ We start with the velocity autocor- 
relation function Z(t), see Eq. |6j that is shown in Fig. 
[T] There are only subtle differences in Z(t) between liq- 
uid and solid phases. For the solid Z(t) is slightly more 
negative for 4 < tui p < 14. However this slight difference 
leads to a much smaller D from the integral in Eq. [7j 
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FIG. 1: (Color on line) Velocity autocorrelation function Z(t) 
versus time t in units of one over the plasma frequency ui p for 
N = 8192 ions at T = 176 for both a liquid configuration (red 
dashed line) and a solid configuration (black solid line). 



Next, Fig. [2] shows histograms of displacements \rj(t+ 
to) — (to) I after a time t = 21000/oj p . These are com- 
puted by simply counting the number of ions that have 
moved a given distance. Figure [2] shows a large peak at 
small distances that corresponds to ions which remain 
near their original lattice site. The width of this peak 
corresponds to thermal oscillations. The amplitude of 
these oscillations are relatively large because the system 
is warm and near the melting temperature. Figure [2] also 
shows smaller peaks at larger distances that correspond 
to ions which have "hopped" one lattice site, two lattice 
sites, etc. Diffusion is seen to be larger for a system that 
started from imperfect crystal initial condition compared 
to a system that started from a perfect bec lattice ini- 
tial condition. We start by presenting additional results 
for perfect body centered cubic lattice initial conditions 
and then we will present results for imperfect crystal and 
amorphous initial conditions. 



How do the ions actually move (hop) from one lattice 
site to the next? This is nontrivial because the system is 
under high pressure and vacancy formation is suppressed. 
Thus there are very few empty sites for the ions to hop 
into. Instead the ions can exchange with their neighbors. 
In Figure [3] we show the final configuration for a small 
3456 ion system that was prepared from perfect bec lat- 
tice initial conditions. Most ions remain near their origi- 
nal lattice site and are shown as small brown dots. These 
ions show oscillations about the lattice sites. However for 
this example, there were 24 ions that moved more than 
1.34a during the finial simulation time of t = 236/w p . 
These ions are shown as larger black disks and are seen 
to be in a ring configuration where ions "hop" to lattice 
sites vacated by other hopping ions. 
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FIG. 2: (Color on line) Histogram of displacements \Tj(t + 
to) — (to) in units of the ion sphere radius a after a time 
t = 21000/cjp. The simulations use N = 8192 ions and are 
at r = 176. The black solid line is the average of 800 con- 
figurations (initial times to) for a system that started from a 
perfect body centered cubic (bec) lattice initial configuration 
while the dashed red line is the average of 8000 configurations 
for a system that started from an imperfect crystal initial con- 
figuration. 



We now present results for the diffusion constant D 
using Eq. [10] with a cutoff parameter R c chosen as the 
location of the minimum in the histograms in Fig. [2] at, 



Rc 



1.07a. 



(12) 



Small changes in this value only lead to slight changes 
in D. To minimize finite size effects we also introduce a 
cutoff range i? cu t in the Yukawa interaction so that Eq. 
[T] becomes 



Vij — ZiZje 2 



-r/\ 



-flcut/A 



Rc. 



e(R c 



(13) 
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FIG. 3: (Color on line) Sample configuration of 3456 ions at 
r = 175. Ions that have moved less than 1.34a in a time 
t — 236/ajp are small brown dots. Ions that have moved more 
than 1.34a are shown as larger black disks and are seen to 
be in a ring configuration where ions "hop" to lattice sites 
vacated by other hopping ions. This system started from a 
perfect bcc lattice. Figure plotted with VMD |46j . 



TABLE I: Diffusion constant D for MD simulations starting 
from perfect body centered cubic lattice initial conditions at 
r = 175. Here D is in units of cj p a 2 with uj p the plasma 
frequency and a the ion sphere radius, N is the number of 
ions, A the screening length, R cu t the cutoff radius in the 
interaction, At the MD time step, t the elapsed time, and 
N con f the number of configurations used to average over the 
inital time to- 
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"This is a continuation of the run described in the line above. 
However, it is at constant energy instead of being at (approxi- 
mately) constant temperature. 



and the potential is zero for r > R cu t- 

We first consider bcc lattice initial conditions. Table U 
presents results for D for different values of A, molecular 
dynamics time step At , number of ions TV, and elapsed 
time t used in Eq. 



10 



We express D in units of uj p 



We find that D increases with decreasing A. For a large 
value of A = 2.70a there are large finite size effects and 
D increases with increasing N. However the increase in 
D in going from N = 27648 to the largest system size 
93312 is small. 

Finite size effects are smaller for the smaller A = 1.82a 
value. Now there is good agreement for N = 27648 and 
93312 and D is only slightly smaller for N = 8192. We 
do not find strong sensitivity to At or t. Furthermore, 
imposing a cutoff on the interaction at large distances 
-Rcut = 8.91 A has only a very small effect on D. For A = 
1.82a and large systems, we find D/uj p a 2 = 1.0 X 10~ 5 . 
As we discuss below, this value, for the solid near the 
melting point T = 175, is about 200 times smaller than 
D for the liquid phase at the same T. 

Diffusion in the solid may involve an energy barrier AE 
since it may be necessary for an ion to pass close to its 
neighbors. This would lead to a temperature dependancc 
D oc Exp(-A£;/T) = Exp(-aT) with d a constant. In 
Table III] we present results for D as a function of T. For 



TABLE II: Diffusion constant D versus V for MD simulations 
using iV = 27648 ions and starting from perfect body centered 
cubic lattice initial conditions. The screening length is A, R cu t 
is the cutoff radius for the interaction, the MD time step is 
59000/0^. 

T A/a Rcut /A D/bj p a 2 



165 1.82 8.91 3.3 x 10~ J 

175 1.82 8.91 1.05 x 10" 5 

185 1.82 8.91 3.3 x 10" 6 

200 1.82 8.91 3.9 x 10" 7 

165 2.70 oo 1.3 x 10" 5 

175 2.70 oo 4.1 x 10" 6 

185 2.70 oo 1.6 x 10~ 6 

200 2.70 oo 1.5 x 10" 7 



r < 185, Table [H*| results are approximately 
fa 6100Exp(-0.115r) 



for A = 1.82a and 
D 



1400 Exp(-0.112T) 



(14) 



(15) 



for A = 2.70a. Note that AE (or d « 0.11) appears to 
be almost independent of A. This would follow if AE is 
dominated by particle interactions at short distances. 
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FIG. 4: (Color on line) Diffusion constant D versus V for 
both the liquid and solid phases. Liquid results are from ref. 
g] while the solid results are for A = 1.82a and assume per- 
fect bcc lattice initial conditions. The dotted red line shows 
metastable (super cooled) liquid results while the blue dashed 
line shows metastable (super heated) solid results. The sys- 
tem is assumed to melt at F = 175. 



In Fig. [4] we plot D as a function of T for both the 
liquid and solid phases. We see that D drops by a large 
factor as the system crystalizes and that D decreases 
much more rapidly, with increasing T, in the solid phase 
compared to the behavior of D in the liquid phase. 

Most of our simulations are at (approximately) con- 
stant temperature where velocities are rescaled every ten 
time steps to keep the kinetic energy fixed. To test the 
sensitivity of our results to this procedure, we have also 
performed a few runs at constant energy, instead of at 
constant temperature, see for example Table [TJ We find 
that D is unchanged within statistics. 



B. Imperfect crystal initial conditions 

We now consider imperfect crystal initial conditions. 
We prepare a liquid initial condition by starting the ions 
off at random positions, with a Maxwell velocity distri- 
bution, and evolving the system at a series of increasing 
r values. The system is observed to equilibrate in a liq- 
uid phase. However as T is increased further the system 
is observed to supercool for T > 175 and then eventually 
freeze. However, often the system freezes into an imper- 
fect crystal with many defects. For example, there can 
be two micro- crystals of different orientations. Once the 
system has frozen, T is decreased back to T = 175 and 
the system is evolved for a long time at this T value and 
the diffusion constant is calculated from Eq. [10] 

Figure [5] shows a configuration of 27648 ions with im- 
perfect crystal initial conditions. Here ions, that have 



moved less than three lattice spacings during the simu- 
lation, are plotted as small gray points while ions, that 
have moved more than three lattice spacings, are plotted 
as large blue spheres. The system froze into two micro- 
crystals of different orientation and the diffusing ions are 
seen to be clustered near the grain boundaries. This sug- 
gests that diffusion in imperfect crystals may be domi- 
nated by motion of the defects rather than by hopping 
of ion chains, such as that shown in Fig. [3j 
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FIG. 5: (Color on line) Configuration of 27648 ions starting 
from imperfect crystal initial conditions. Ions that move only 
a small distance are small gray points. Ions that have moved 
over three lattice spacings, during the simulation time of t = 
59000/cjp, are shown as large blue spheres. These are seen to 
be clustered at the grain boundaries. The initial conditions 
included two micro-crystals of different orientation. Figure 
plotted using VMD gB]. 



In Table III we present results for D for imperfect crys- 
tal initial conditions. Note that lines 3 and 4 in Table Hill 
and lines 5 and 6 correspond to independently prepared 
initial conditions. There is some variation in results for 
different simulations. This may reflect differences in the 
number and kind of defects present in the initial condi- 
tions. We see that D for imperfect crystal initial condi- 
tions is two to four times larger than D for perfect bcc 
lattice initial conditions. We also see that D may be less 
sensitive to the screening length for imperfect crystal ini- 
tial conditions. 

It is possible that D will evolve slowly with simula- 
tion time to for these imperfect crystal simulations. Note 
that we do not find rapid variation of D with t . How- 
ever, we have not attempted to determine how D might 
evolve over long times by continuing an imperfect crystal 
simulation for very long times. 

The final simulation listed in Table |III| was prepared 
by very slowly cooling a liquid configuration that started 
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TABLE III: Diffusion constant D for MD simulations starting 
from imperfect crystal initial conditions at T — 175. Here D 
is in units of ui p a 2 with uj p the plasma frequency and a the ion 
sphere radius, N is the number of ions, A the screening length, 
Rent the cutoff radius in the interaction, At the MD time step, 
t the elapsed time, and iV COI1 fig the number of configurations 
used to average over the initial time to. 



N 


A/a 


^cut 


/A Atui p 


tU) p 


-^V config 


D/ujpO, 2 




3456 


2.70 


OO 


0.047 


4700 


380 


4.9 x 10" 


5 


8192 


2.70 


OO 


0.12 


9600 


8000 


1.4 x 10" 


-5 


27648 


2.70 


OO 


0.12 


59000 


400 


8.9 x 10" 


-6 


27648° 2.70 


OO 


0.12 


59000 


400 


8.8 x 10" 


-6 


27648 


2.70 


OO 


0.12 


59000 


800 


1.9 x 10" 


-5 


27648 


1.82 


8.91 


0.12 


59000 


500 


2.3 x 10" 


-5 


27648 


1.82 


OO 


0.06 


472000 401 


1.1 x 10" 


-5 



"This is a continuation of the run described in the line above. 
However it is at constant energy instead of being at (approximately) 
constant temperature. 



TABLE V: Approximate time to for amorphous systems 
to crystalize, after the systems have been instantaneously 
quenched from F = 175 to different F values, see text. The 
number of ions is N. 

n r t up 

8192 500 24,000 
8192 600 47,000 
8192 700 142,000 
8192 1000 240,000 
8192 1500 390,000 
8192 1750 > 25,000,000 
27648 500 400,000 
27648 1000 > 6, 000, 000 

some variation in these results for D depending on the 
number and kind of defects present in the initial condi- 
tions. 



TABLE IV: Diffusion constant D versus F for MD simulations 
using N — 27648 ions and starting from a single imperfect 
crystal initial condition. The screening length is A, R cu t is 
the cutoff radius for the interaction, the MD time step is 
Atcup = 0.12 and t = 59000a> p . 

T A/a R cut /X D/ujpa 2 

175 1.82 8.91 2.3 x 10" 5 

185 1.82 8.91 7.2 x 10" 6 

200 1.82 8.91 2.1 x 10" 6 

225 1.82 8.91 3.5 x 10" 7 

250 1.82 8.91 5.8 x 10" 7 



at T = 175, at a rate of dT/dt = 2.1 x lfr 4 u;p, until the 
configuration froze at T = 283. The resulting solid config- 
uration was then heated back up to T = 175. Finally, the 
system was evolved at T = 175 for a time to = 59000/cj p 
before taking D data. This system was observed to be a 
nearly perfect bec lattice, and the value for D in Table 
|III| agrees with our results for nearly perfect bec lattices 
in Section III A This strongly suggests that white dwarf 
and neutron star plasmas will freeze into nearly perfect 
bec crystals, since any astrophysical cooling time scale is 
likely very much longer than this MD cooling time scale. 
This is also consistent with our results in Section MI CI 
for amorphous systems, see below. 

In Table |IV| we present results for D versus T for a 
single imperfect crystal initial condition. Here D was 
calculated at T — 175. Next the velocities of the final 
r = 175 configuration were rescaled to T = 185 and the 
system was equilibrate at T = 185 and D determined. 
This process was repeated for larger T. We see that D 
decreases with increasing T far more slowly than does D 
for a bec lattice. This suggests that D is dominated by 
the motion of defects and that these defects continue to 
move even at low temperatures where the ion hopping 
shown in Fig. [3] is very unlikely. Note that we expect 



C. Amorphous initial conditions 

In this subsection, we present results for amorphous 
initial conditions. The imperfect crystal initial condi- 
tions in Subsection IIIIBI involved a small amount of su- 
percooling, until a configuration froze. We now consider 
much greater supercooling. We start with a liquid config- 
uration of N=8192 ions, that is equilibrated at T = 175. 
The screening length is A = 1.82a, R cu t = 8.91A, and the 
time step is Atuj p = 0.12. We quench the system instan- 
taneously to a large T value by rescaling the velocities, 
and then we evolve the resulting amorphous system at 
(approximately) constant temperature until the system 
largely crystalizes. Note that quenched initial configu- 
rations for different T values were prepared by rescaling 
the velocities of the same T = 175 liquid configuration. 
Table |V| lists the time needed to crystalize for different T 
values. This time increases with increasing T (amount of 
supercooling). However, this time only increases approx- 
imately linearly with T for V < 1500. This suggests that 
diffusion is relatively fast in the amorphous system, and 
that the amorphous to crystal transition does not involve 
a large energy barrier. We find that the system is able to 
crystalize, even at T = 1500 where the temperature is 8.6 
times lower than the melting temperature. However a fi- 
nal run that was quenched to T = 1750 was not observed 
to crystalize before a time 25, 000, OOO/wp. We refer to 
these quenched systems as amorphous. However, it may 
be more appropriate to call them polycrystalline because 
they are observed to have many small crystal domains 
of different orientation. These polycrystalline states are 
observed to undergo rapid transitions to single crystals 
except at the largest F, see below. 

To study finite size effects we now consider larger sys- 
tems with TV = 27648 ions. We start with a liquid config- 
uration equilibrated at T = 175, and quench the system 
instantaneously to T = 500. Figure [6] shows D versus 



simulation time to- The diffusion constant D first de- 
creases rapidly with time as the quenched system an- 
neals. Next D remains more or less constant for a long 
time. Suddenly near to — 400, 000/cj p there is a large 
spike in D. We calculate D with both t — 59000/w p and 
4720/wp. The larger t gives D with less statistical noise, 
while the smaller t gives better time resolution and shows 
that the event near to — 400, OOO/wp is very rapid. 




FIG. 6: (Color on line) Diffusion constant D versus simulation 
time to, see Eq. |10| for an amorphous system of N = 27648 
ions at V — 500. The diffusion constant D is calculated using 
a time difference of t = 59, 000 /ajp (solid black line) and t = 



10 The sample was prepared 



/ui p (dotted red line) in Eq. 
by instantaneously quenching a liquid from T = 175 to T = 
500. 



The configuration of the system just before the event is 
shown in Fig. [7] The system is seen to be in a polycrys- 
talline state with many small crystal domains. Figure 
[8] shows the configuration of the system just after the 
event. Now the system is an imperfect single crystal. 
This demonstrates that diffusion is fast enough, at least 
at r = 500, for the system to crystalize. Finally in Fig. 
[9j we show D as a function of simulation time to for an 
TV = 27648 ion system quenched to T = 1000. Again D 
starts off large and decreases rapidly as the system starts 
to equilibrate. Three large peaks are observed in D near 
t = 150, 000, 400, 000, and 5x 10 6 /^ p . These correspond 
to events where small micro-crystals rearrange and grow 
and the bond angle metric Qe increases, as we discuss 
below, see Fig. [TUJ However the system is still polycrys- 
talline after the events. 

To quantify the crystalline order in these simulations, 
we consider a metric based on bond angles 47-49 , see 
also [30] . Ion i is said to be bonded to ion j if it is within 
a distance b — 2.44a that corresponds to a minimum 
in the radial distribution function g{r). This distance is 
chosen to include the eight nearest neighbors and six next 
nearest neighbors in a perfect body centered cubic lattice. 
Let Oij and <f>ij be the polar and azimuthal angles of the 



radius from ion i to ion j. 
harmonic, 



We calculate the spherical 



'At. 



lm \yij i ) i 



(16) 



and average over all ~ 14./V bounds for a given configu- 
ration, 



Qlm = (Qlm(fij)} ■ 



(17) 



This quantity depends on the orientation of a crystal lat- 
tice with respect to the simulation volume. Therefore, we 
calculate the rotationally invariant quantity Qi [13 2PJ , 



Qi = 



4tt 



2/ 



i 

m— — l 



1/2 



(18) 



This provides a measure of the crystalline order of a con- 
figuration. In general, Qi is small for a liquid or amor- 
phous configuration and Qi is large for a perfect crystal. 
Our calculations of Qi, for a range of even I, show that 
Qe is most sensitive to crystalline order. We find that 



Qe = 0.51069 



(19) 



for a perfect bcc lattice and Qe = 0.57452 for a per- 
fect face centered cubic lattice, see also [H]. Note that 
Qi is small for odd I. In Fig. [10] we show Qe versus 
simulation time to- In general, Qe grows with to- How- 
ever, the amount of time necessary for Qe to grow can 
increase strongly with increasing system size N or T. A 
plateau near Qe ~ 0.17 is seen for all four systems in 
Fig. 10 This suggests a possible metastable intermedi- 
ate state. The simulations with N = 27648 at T = 500 
and iV = 8192 at T = 1500 show a rapid rise in Q 6 , near 
to = 4 x 10 5 /a;p, during transitions to single crystals. 
For N = 27648 at T = 1000 and N = 8192 at T = 1750, 
Qe is increasing with time. These systems have not yet 
evolved to single crystals. However, the continued rise of 
Qe with time strongly suggests that these systems will 
evolve to single crystals at later times. In summary, the 
continued rise of Qe with time, as shown in Fig. |10| 
demonstrates that these quenched systems are evolving 
with time towards single crystals, and that they are un- 
likely to remain amorphous. 

We conclude that an amorphous solid will not form 
even with large amounts of super cooling, where the tem- 
perature is rapidly quenched by up to a factor of 10 be- 
low the melting temperature. Instead, diffusion is fast 
enough so that the system will form a regular crystal. 
Our results strongly suggest that Coulomb solids in the 
interior of cold white dwarf stars and the crust of neu- 
tron stars will be crystalline and not amorphous. This 
is consistent with observations of rapid crust cooling of 
neutron stars following extended periods of accretion [2"4T - 
[27] . This rapid cooling implies a high crust thermal con- 
ductivity, that agrees with the conductivity of a regular 
crystal, and is larger than the conductivity expected for 
an amorphous solid. 
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FIG. 7: (Color on line) Configuration of 27648 ions at F = 500 
after a simulation time to = 350, 000/cj p . The sample was 
prepared by instantaneously quenching a liquid from F = 175 
to T = 500. Figure plotted using VMD gg]. 




FIG. 8: (Color on line) Configuration of 27648 ions at F = 500 
after a simulation time to = 450, 000/cj p . The sample was 
prepared by instantaneously quenching a liquid from F = 175 
to T = 500. Figure plotted using VMD [46]. 



IV. CONCLUSIONS 

Diffusion in coulomb crystals can be important for the 
structure of the crust of neutron stars. In this paper, 
we perform molecular dynamics simulations of one com- 
ponent coulomb crystals to study the diffusion constant 
D. We find that D is non-zero, at least near the melt- 
ing temperature, and that D for Coulomb crystals with 
relatively soft-core 1/r interactions is in general larger 




LCD 

p 

FIG. 9: (Color on line) Diffusion constant D versus simu- 
lation time to for a 27648 ion system at V = 1000, using 
t = 59000/cjp. The vertical lines mark diffusion features that 
are also indicated in Fig. |10| The system was prepared by 
instantaneously quenching a liquid from T = 175 to F = 1000. 



than D for Lennard- Jones or other solids with harder- 
core (more singular) interactions. 

We find that diffusion, for simulations that start from 
a perfect body-centered-cubic lattice, involves the ex- 
change of ions in ring-like configurations. Here ions 
"hop" in unison with one ion replacing another with- 
out the formation of long lived vacancies. This may be 
true because vacancy formation is strongly suppressed 
because of the large pressure. The diffusion constant 
D decreases rapidly, for temperatures below the melt- 
ing point, suggesting that these ring-like configurations 
have a high activation energy. 

We also calculate diffusion for simulations that start 
from imperfect crystal initial conditions. Here a liq- 
uid configuration, at a temperature somewhat below the 
melting point, spontaneously freezes to an (in general) 
imperfect crystal that may contain defects such as dis- 
locations and grain boundaries. Note that these con- 
figurations involve one or more micro-crystals and are 
not amorphous. For these systems, D is larger than D 
for perfect bcc lattice configurations and decreases more 
slowly with decreasing temperature. This suggests that 
D for imperfect crystals is dominated by the motion of 
the crystal defects rather than the hopping of ions in a 
perfect crystal. Therefore, observations of D may help 
characterize the imperfections in a Coulomb crystal. 

Finally, we studied diffusion in "amorphous" systems 
where the temperature was instantaneously quenched to 
much lower values. We find that D is large. Indeed most 
of our amorphous simulations are observed to sponta- 
neously transform to either a single crystal or a small 
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': r=1500, N=8192 


T=500, N=27648 




r=1000, N=27648 ; 




aarT^n^u 
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FIG. 10: (Color on line) Bond angle metric Qe, see Eqs. 
|16|17|18"1 versus simulation time to for amorphous systems 
that were instantaneously quenched from a F = 175 liquid 
at to = 0. Fhe number of ions in the simulation N and 
coulomb parameter F are indicated. Fhe vertical red lines, 
for r = 1000 and N = 27648, indicate diffusion features that 
are seen in Fig. [9] 



number of crystal domains. This strongly suggests that 
Coulomb solids in white dwarf and neutron stars are crys- 
talline, rather than amorphous. This is in agreement with 
X-ray observations of rapid neutron star crust cooling 
that imply a large thermal conductivity. 

It is an important open problem to determine the equi- 
librium distribution of defects in a Coulomb crystal. It 
may be difficult to determine this directly from molecu- 
lar dynamics simulations because it can take a very long 
time for defects to equilibrate. However, we find that 
diffusion is relatively fast. This suggests that astrophysi- 
cal coulomb solids will have had plenty of time to anneal 
to nearly perfect crystals with relatively few defects. Fi- 
nally, the diffusion constant that we find for a pure bec 
lattice may provide a lower limit on D for an equilibrated 
system. This is because the presence of defects is only 
expected to increase D over that for a perfect crystal. 
In future work we plan to study D for multicomponent 
Coulomb solids. For a given species i, we expect a rich 
behavior for the diffusion constant Di depending on how 
the charge of an ion Zj compares to the average charge 
of the ions that make up the crystal lattice. 

We thank Andrey Chugunov for very helpful com- 
ments. This research was supported in part by DOE 
grant DE-FG02-87ER40365 and by the National Science 
Foundation through TeraGrid resources provided by Na- 
tional Institute for Computational Sciences under grant 
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